09 - Remote Sensing with Landsat package

Working with satellite imagery in R

In this exercise you learn to use the functions of the landsat package to process satellite imagery, specifically, how to:

  • streamline histograms so you can combine images
  • manipulate multi-band imagery,
  • extract and create new data from images, such as NDVI and SAVI
  • classify image values using kmeans unsupervised algorithm to detect similar areas
  • segment features in satellite imagery

Task 1: Set up your workspace

Start by installing the required packages first. rasterVis and RColoBrewer are used for visualisation of rasters, lattice and latticeExtra for extra graphical utilities, while landsat provides the imagery samples and rgl allows for 3D visualisation.

# install.packages(c('rasterVis', 'RColorBrewer', 'landsat', 'lattice',
# 'latticeExtra', 'rgl', 'itcSegment'))

library(lattice)
library(latticeExtra)
library(RColorBrewer)
library(rasterVis)
# library(rgdal)
library(rgl)

Task 2: Pre-processing of Landsat datasets

Landsat packages offers sample Landsat satellite imagery decomposed into single bands, and labelled by the month of capture (25 November 2002). You will load and practice raster analysis with these 300x300 pixel samples. Note that each single band image shows as panchromatic (greyscale) with values ranging from 0 (white) to 255 (black), with the interim colors being shades of grey. This is unsigned 8 bit imagery. The number behind the filename (e.g. nov3) refers to the color band of the image: 1 = Green, 2 = Blue, 3= Red, 4 = Near-infrared. Once you combine these color bands and plot with plotRGB(), you can see the true- or false-color imagery depending on which color you load into which RGB channel.

The images are in the SpatialGridDataFrame format and need to be converted to raster format before manipulation.

library(landsat)

## load indvidual band image data from landsat package
`?`(nov)

# load band#3 red channel of the image
data(nov3)
plot(nov3)

data(nov4)
plot(nov4)

Task 3: Load elevation data and plot it in 3D

dem in the landsat package denotes a ‘digital elevation model’. Once you load it you can convert it to Formal Raster Layer with raster() function and then continue processing using the usual raster package functions.

plot3D() is a neat function in the rasterVis package, which opens a separate windows and plots elevation values in 3D space if present in the raster. You need to close the window if you want to update or plot another object.

# Load and plot the digital elevation model in landsat package
library(raster)
data(dem)
plot(dem)

dem <- raster(dem)

Plot in 3D

# RasterVis package function plot3D() loads a neat 3D viewer, which opens in a
# new window If you don't see the plot3D() function working, enable the
# rglwidget() by uncommenting the next line

# rglwidget()
plot3D(dem, rev = T, zfac = 1)

Task 4: Load and explore RGB data components

# let's load data for July and explore it
data("july1")
data("july2")
data("july3")
data("july4")
j1 <- raster(july1)  # blue
j2 <- raster(july2)  # green
j3 <- raster(july3)  # red
j4 <- raster(july4)  # near-infrared

## check out the image histogram\t
plot(j1)

plot(j4)

hist(j1, main = "Band 1 - Blue of Landsat")

boxplot(j1, main = "Band 1 - Blue of Landsat")

Task 5: Plot RGB image

# take the June data and create an RGB image and drape it over the 3D model

### Reorder to R - G - B and create a multi-layer rasterBrick and also a
### false-color rasterStack!
myRGB <- brick(j3, j2, j1)  # brick creates new object
myCIR <- stack(j4, j3, j2)  # stack stores connections only


### let's see how the NIR, R, and G bands relate (from lattice)
splom(myCIR, varname.cex = 1)  # scatter plot matrix

# let´s plot in full colour!
plot(myRGB)

## better
plotRGB(myRGB)

plotRGB(myCIR)

Ok, you can probably see something?

Task 6: Manipulate image rendering by histogram stretching

First, let’s use histogram stretch

## different stretches - here histogram based
plotRGB(myRGB, stretch = "hist")

plotRGB(myCIR, stretch = "hist")

Next, let’s try linear stretch

## different stretches - here linear stretch
plotRGB(myRGB, stretch = "lin")

plotRGB(myCIR, stretch = "lin")  # in CIR red = green!

Any idea what the red color represents in myCIR?

Finally, drape one of the images over a 3D model. This bit can be a bit touchy and take a while to get to work. I tend to run the plot3D() lines alone to generate the 3D view in a pop-up window rather than rmarkdown output. Everytime you wish to refresh the view, you must close the pop-up window.

## finally...
rglwidget()  # this widget helps get the first view rendered in rmarkdown, refreshing however is more tedious

`?`(rglwidget())

t <- plot3D(dem, col = rainbow(255))  ## you need to close RGL device manually first and then run this line!


# Save it to a file.  This requires pandoc
filename <- tempfile(fileext = ".html")
htmlwidgets::saveWidget(rglwidget(), filename)
browseURL(filename)

plot3D(dem, drape = j4)  ## should drape image j4 over DEM, if problematic try in .R script and watch for a pop-up widget. Ask Adela to demo!

Additionally, the histograms might need some adjustment to balance/equalize them before plotting. Notice how the November and July value scales are different at 80 and 250 max respectively. Clearly there is much more variation in July growth and we need to stretch the November histogram a bit to have a comparable view.

# histogram/color adjustments
data(nov3)
data(july3)
par(mfrow = c(2, 2))
image(nov3, main = "nov")
image(july3, main = "july")
plot(nov3, main = "nov")

plot(july3, main = "july")

# let's create a matching histogram by stretching nov3 to july3
nov3.newH <- histmatch(master = july3, tofix = nov3)
image(nov3.newH$newimage, main = "new nov")  # look at it on the fly

# write the new image with a new histogram
n3new <- raster(nov3.newH$newimage)

# convert existing images to raster format before plotting
n3 <- raster(nov3)
j3 <- raster(july3)

# plot (the mfrow is set up to work in R script, so just click through the
# rmarkdown view)
par(mfrow = c(1, 3))

hist(n3, main = "Nov")
hist(j3, main = "July")
hist(n3new, main = "new Nov")

#### most important corrections are atmospheric and topographic however, these
#### are too complex to cover here...see package help

par(mfrow = c(1, 1))  # remember to change plotting to a single window if in .R script

# Plot the equalized-histogram images
plot(raster(july3), main = "july")

plot(n3new, main = "new november")

Well done on equalizing histograms over the same area. Now you can see different kinds of vegetation thriving at different times of year.

Task 7: Create new information from satellite imagery: NDVI

Let us calculate the Normalized Difference Vegetation Index (NDVI) and see where the vegetation grows most in our Landsat image: Remember the formula for NDVI is: (NIR - RED) / (NIR + RED)

# prep imagery
n3 <- raster(nov3)  # RED
n4 <- raster(nov4)  # NIR

ndvi <- (n4 - n3)/(n4 + n3)
ndvi
class      : RasterLayer 
dimensions : 300, 300, 90000  (nrow, ncol, ncell)
resolution : 30, 30  (x, y)
extent     : 390045, 399045, 4482105, 4491105  (xmin, xmax, ymin, ymax)
crs        : NA 
source     : memory
names      : layer 
values     : -0.3114754, 0.5664336  (min, max)
plot(ndvi)

# uncomment and run in 3D
plot3D(dem, drape = ndvi, zfac = 1.5)

### remove values below zero
ndvi[ndvi <= 0] <- NA
plot(ndvi)

## plot again plot3D(dem, drape=ndvi)

## different way to plot in 2D
library(rasterVis)
levelplot(ndvi, col.regions = bpy.colors(100))

Now you can see the areas of the highest reflectance and thus most healthy vegetation in November

—- skip to task 9 if short on time

Task 8: Create new information from satellite imagery: SAVI

SAVI stands for Soils Adjusted Vegetation Index, and this is another calibrated view of the ground.

### another index SAVI (soil adjusted vegetation index)
ndvi <- (n4 - n3)/(n4 + n3)
savi <- (n4 - n3)/((n4 + n3) * 0.25)  # with L=1 -> similar to NDVI

### let´s compare visually
par(mfrow = c(1, 2))
plot(savi, main = "SAVI")
plot(ndvi, main = "NDVI")

par(mfrow = c(1, 1))

Task 9: Unsupervised Classification with k-means

We would like to isolate and better see the clusters of growth within our image. We will run kmeans function on a composite image in order to cluster data based on similarity or similar groups!

# first, let's select an image and make it into a brick including ndvi
data(nov2)
data(nov1)
n2 <- raster(nov2)
n1 <- raster(nov1)
ndvi
class      : RasterLayer 
dimensions : 300, 300, 90000  (nrow, ncol, ncell)
resolution : 30, 30  (x, y)
extent     : 390045, 399045, 4482105, 4491105  (xmin, xmax, ymin, ymax)
crs        : NA 
source     : memory
names      : layer 
values     : -0.3114754, 0.5664336  (min, max)
which(is.na(as.data.frame(ndvi)))  # should be 0, rerun the ndvi creation if not.
integer(0)
# create a new composite brick out of the available data
myNewBrick <- brick(n4, n3, n2, n1, ndvi)
splom(myNewBrick)

plot(myNewBrick)

Next, run the kmeans classification. Beware that the kmeans function does not tolerate NA/INF/NaN and similar values. Our new brick should not have any but in future classification remember that you need to get around them, either by exclusion or substitution via mean values.

# Run kmeans classification on the values in your new brick Read on
# Thresholding here:
# https://rspatial.org/raster/rs/3-basicmath.html#vegetation-indices
ICE_df <- as.data.frame(myNewBrick)
set.seed(99)


cluster_ICE <- kmeans(ICE_df, 10)  ### kmeans, with 4 clusters
str(cluster_ICE)
List of 9
 $ cluster     : int [1:90000] 4 4 4 6 6 7 6 6 4 1 ...
 $ centers     : num [1:10, 1:5] 77.8 33.1 93.3 66.8 44.9 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:10] "1" "2" "3" "4" ...
  .. ..$ : chr [1:5] "X.data.LANDSAT.MASTER.Gorskiout.G20021125_7_band4.geo.asc" "X.data.LANDSAT.MASTER.Gorskiout.G20021125_7_band3.geo.asc" "X.data.LANDSAT.MASTER.Gorskiout.G20021125_7_band2.geo.asc" "X.data.LANDSAT.MASTER.Gorskiout.G20021125_7_band1.geo.asc" ...
 $ totss       : num 20611557
 $ withinss    : num [1:10] 169596 158034 131141 241482 231449 ...
 $ tot.withinss: num 1874040
 $ betweenss   : num 18737517
 $ size        : int [1:10] 4612 10833 1783 6804 16124 8862 16155 3728 4620 16479
 $ iter        : int 10
 $ ifault      : int 0
 - attr(*, "class")= chr "kmeans"
# cluster_ICE <- cluster::clara(ICE_df, 4) ### another option, clara, with 4
# clusters

# convert cluster information into a raster for plotting
clusters <- raster(myNewBrick)  ## create an empty raster with same extent than ICE
clusters <- setValues(clusters, cluster_ICE$cluster)  # convert cluster values into raster
clusters
class      : RasterLayer 
dimensions : 300, 300, 90000  (nrow, ncol, ncell)
resolution : 30, 30  (x, y)
extent     : 390045, 399045, 4482105, 4491105  (xmin, xmax, ymin, ymax)
crs        : NA 
source     : memory
names      : layer 
values     : 1, 10  (min, max)
plot(clusters)

# uncomment to plot the clusters in 3D over the DEM
plot3D(dem, drape = clusters, col = c("white", "green", "blue", "yellow"))

# calculate the average spectral signature of 1-4 bands of growth
ICE_mean <- zonal(myNewBrick, clusters, fun = "mean")
ICE_mean  # see the values for ndvi (layer) being most distinct
      zone X.data.LANDSAT.MASTER.Gorskiout.G20021125_7_band4.geo.asc
 [1,]    1                                                  77.78751
 [2,]    2                                                  33.07745
 [3,]    3                                                  93.29893
 [4,]    4                                                  66.75720
 [5,]    5                                                  44.89792
 [6,]    6                                                  57.19567
 [7,]    7                                                  50.56335
 [8,]    8                                                  56.26502
 [9,]    9                                                  47.30498
[10,]   10                                                  39.66339
      X.data.LANDSAT.MASTER.Gorskiout.G20021125_7_band3.geo.asc
 [1,]                                                  41.72029
 [2,]                                                  31.33093
 [3,]                                                  40.36960
 [4,]                                                  43.96326
 [5,]                                                  37.63799
 [6,]                                                  42.74362
 [7,]                                                  40.45453
 [8,]                                                  50.90504
 [9,]                                                  45.32706
[10,]                                                  34.33989
      X.data.LANDSAT.MASTER.Gorskiout.G20021125_7_band2.geo.asc
 [1,]                                                  45.40590
 [2,]                                                  34.81261
 [3,]                                                  46.12900
 [4,]                                                  45.14065
 [5,]                                                  38.55848
 [6,]                                                  42.28154
 [7,]                                                  39.87242
 [8,]                                                  48.89753
 [9,]                                                  44.43831
[10,]                                                  36.50592
      X.data.LANDSAT.MASTER.Gorskiout.G20021125_7_band1.geo.asc      layer
 [1,]                                                  58.13964 0.30225062
 [2,]                                                  52.30970 0.02607170
 [3,]                                                  58.31408 0.39513993
 [4,]                                                  58.60553 0.20678500
 [5,]                                                  54.86039 0.08794067
 [6,]                                                  56.87599 0.14495138
 [7,]                                                  55.61634 0.11117963
 [8,]                                                  62.34576 0.04988818
 [9,]                                                  59.38939 0.02106308
[10,]                                                  53.31756 0.07190780

Note that you have aggregated the final cluster raster by using the focal() function using the mean function on the four clusters identified by kmeans() as similar.

Task 10: What is the trend in de-/afforestation? - Individual tree crown segmentation

The ITC (Individual Tree Crowns) delineation approach finds local maxima within imagery that contains subtle color differences, such as the canopy image provided. The itcIMG() function designates these maxima as tree tops, then uses a decision tree method to grow individual crowns around the local maxima.

The image we use is based on LiDAR (Light Detection and Ranging) in xyz format. To get the segmentation going you will need either rgeos with itcSegment version 07 or the terra and current itcSegment packages

# install.packages('terra') install.packages('itcSegment')

library(terra)
library(itcSegment)

imgData <- rast("../data/imgData.tif")

plot(imgData)

# Use the itcIMG() function to detect and grow the individual crowns

se <- itcIMG(imgData, epsg = 32632)
summary(se)
       ID               X                Y               CA_m2      
 Min.   :  1.00   Min.   :676744   Min.   :5091658   Min.   : 6.48  
 1st Qu.: 30.25   1st Qu.:676754   1st Qu.:5091661   1st Qu.:34.12  
 Median : 59.50   Median :676776   Median :5091663   Median :43.74  
 Mean   : 59.50   Mean   :676782   Mean   :5091663   Mean   :45.66  
 3rd Qu.: 88.75   3rd Qu.:676802   3rd Qu.:5091666   3rd Qu.:56.17  
 Max.   :118.00   Max.   :676832   Max.   :5091669   Max.   :95.58  
# What is the product of the function? Is it a raster or vector?
se
 class       : SpatVector 
 geometry    : polygons 
 dimensions  : 118, 4  (geometries, attributes)
 extent      : 676744.4, 676837.1, 5091659, 5091747  (xmin, xmax, ymin, ymax)
 coord. ref. : WGS 84 / UTM zone 32N (EPSG:32632) 
 names       :    ID         X         Y CA_m2
 type        : <int>     <num>     <num> <num>
 values      :     1 6.768e+05 5.092e+06 59.94
                   2 6.768e+05 5.092e+06 82.62
                   3 6.768e+05 5.092e+06 54.27
plot(se, axes = T)

### Let´s overlay the image and the product of segmentation (run both lines)
plot(imgData)
plot(se, axes = F, add = T)

Task 11: Visualise the segmentation result in Leaflet

You can probably do all of this yourself, but here is a hint about projecting the SpatialPolygonsDataFrame, just in case:

# What are we reprojecting and what to? geographical coordinates or?
se  # it is a Spatial object
 class       : SpatVector 
 geometry    : polygons 
 dimensions  : 118, 4  (geometries, attributes)
 extent      : 676744.4, 676837.1, 5091659, 5091747  (xmin, xmax, ymin, ymax)
 coord. ref. : WGS 84 / UTM zone 32N (EPSG:32632) 
 names       :    ID         X         Y CA_m2
 type        : <int>     <num>     <num> <num>
 values      :     1 6.768e+05 5.092e+06 59.94
                   2 6.768e+05 5.092e+06 82.62
                   3 6.768e+05 5.092e+06 54.27
crs(se)  # what is its crs?
[1] "PROJCRS[\"WGS 84 / UTM zone 32N\",\n    BASEGEOGCRS[\"WGS 84\",\n        DATUM[\"World Geodetic System 1984\",\n            ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n                LENGTHUNIT[\"metre\",1]]],\n        PRIMEM[\"Greenwich\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        ID[\"EPSG\",4326]],\n    CONVERSION[\"UTM zone 32N\",\n        METHOD[\"Transverse Mercator\",\n            ID[\"EPSG\",9807]],\n        PARAMETER[\"Latitude of natural origin\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8801]],\n        PARAMETER[\"Longitude of natural origin\",9,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8802]],\n        PARAMETER[\"Scale factor at natural origin\",0.9996,\n            SCALEUNIT[\"unity\",1],\n            ID[\"EPSG\",8805]],\n        PARAMETER[\"False easting\",500000,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8806]],\n        PARAMETER[\"False northing\",0,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8807]]],\n    CS[Cartesian,2],\n        AXIS[\"(E)\",east,\n            ORDER[1],\n            LENGTHUNIT[\"metre\",1]],\n        AXIS[\"(N)\",north,\n            ORDER[2],\n            LENGTHUNIT[\"metre\",1]],\n    USAGE[\n        SCOPE[\"Navigation and medium accuracy spatial referencing.\"],\n        AREA[\"Between 6°E and 12°E, northern hemisphere between equator and 84°N, onshore and offshore. Algeria. Austria. Cameroon. Denmark. Equatorial Guinea. France. Gabon. Germany. Italy. Libya. Liechtenstein. Monaco. Netherlands. Niger. Nigeria. Norway. Sao Tome and Principe. Svalbard. Sweden. Switzerland. Tunisia. Vatican City State.\"],\n        BBOX[0,6,84,12]],\n    ID[\"EPSG\",32632]]"
# Project the SpatialPolygon using the spTransform() function se4326 <-
# spTransform(se,CRS( '+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
# '))

# Project the SpatialPolygon using sf library
library(sf)
se4326 <- st_transform(st_as_sf(se), 4326)

# ...then we can combine them with leaflet
library(leaflet)
leaflet() %>%
    addTiles() %>%
    addRasterImage(imgData) %>%
    addPolygons(data = se4326, weight = 1, color = "black")  # Tadaa

Can you figure out where this forested landscape is from?

Task 12: OPTIONAL: a more demanding segmentation example with a bigger raster!

And because we have it in the data/ folder, try with this larger example.

r <- rast("../data/myDem_subset.tif")
r
plot(r)

r.se <- itcIMG(r, epsg = 25829, ischm = T)  # can be a bit slow (2-3mins)
summary(r.se)
plot(r)
plot(r.se, axes = F, add = TRUE)

# Adjust 'th' argument for excessive capture of small growth
r.se5 <- itcIMG(r, epsg = 25829, th = 5, ischm = T)  # th - how low should algorithm be looking for canopy
plot(r)
plot(r.se5, axes = T, add = TRUE)

# Write the result to shapefile
td
st_write(st_as_sf(r.se), "../data/itcTrees_subset", driver = "ESRI Shapefile")

# want to see it in Leaflet?
library(sf)
rse <- st_read("../data/itcTrees_subset/itcTrees_subset.shp")
rse4326 <- st_transform(rse, crs = 4326)

# Control question: where is this landscape from?
leaflet() %>%
    addTiles() %>%
    addProviderTiles("Esri.WorldPhysical") %>%
    # addProviderTiles('Esri.WorldImagery') %>%
addRasterImage(r) %>%
    addPolygons(data = rse4326, weight = 1, color = "black")  # Neat :)

The End

Similar approaches can be used when mapping socio-cultural phenomena in satellite imagery, from mosaicing images of night lights as proxies of economic performance, or detecting phenomena in the landscape such as burial mounds, growing urban sprawl, or tracing the outlines of scanned line drawings. (In the latter two you may need to base the classification on reflectance or edge detection rather than elevation.)